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Abstract. In the second magnctohydrodynamic (MHD) ballooning stable 
domain of a high-beta tokamak plasma, the Schrodinger equation for ideal MHD 
shear Alfven waves has discrete solutions corresponding to standing waves trapped 
between pressurc-gradicnt-induced potential wells. Our goal is to understand how 
these so-called a-induced toroidal Alfven eigenmodes (aTAE) are modified by the 
effects of finite Larmor radii (FLR) and kinetic compression of thermal ions in the 
limit of massless electrons. In the present paper, we neglect kinetic compression 
in order to isolate and examine in detail the effect of FLR terms. After a review of 
the physics of ideal MHD oTAE, the effect of FLR on the Schrodinger potential, 
eigenfunctions and eigenvalues are described with the use of parameter scans. 
The results are used in a companion paper to identify instabilities driven by 
wave-particle resonances in the second stable domain. 



1. Introduction 

The second magnctohydrodynamic (MHD) ballooning stable domain [? is an 
attractive parameter regime which may be utilized to achieve high plasma pressures 
in toroidal magnetic confinement devices such as tokamaks with the goal to create 
thermonuclear fusion conditions. Experimental access to the second stable domain 
remains a challenging task requiring reliable profile control (to circumvent the 
ballooning unstable domain via negative magnetic shear and the average magnetic 
well), optimized shaping of the plasma cross-section (to maximize the average magnetic 
well), and utilization of potentially stabilizing mechanisms (e.g., sheared flows [? ]). 

In the meantime, theoretical and numerical studies have advanced to explore the 
properties of the second stable domain. The present paper and its companion [? ] are 
part of this effort and were motivated by a discovery made by Hu & Chen [? ] : using 
the s-a model equilibrium [? ], where s is the magnetic shear and a the normalized 
pressure gradient, Hu & Chen have shown that the second ballooning stable domain 
is populated by discrete ideal MHD Alfven eigenmodes. These so-called a-induced 
toroidal Alfven eigenmodes (aTAE) are standing waves trapped between potential 
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barriers. These barriers represent the modulation of magnetic curvature and shear 
due to the effect of a steep pressure gradient (corresponding to large a) and cause 
reflections of shear Alfven waves (SAW). 

Hu & Chen [? ] have further demonstrated that aTAEs exist for negative 
magnetic shear and that these modes may be resonantly excited through interactions 
with a sparse energetic ion population. These results arc based on a hybrid model 
where an ideal MHD bulk plasma interacts with gyrokinctic energetic ions. Our goal 
is to extend this earlier study by including non-ideal effects of a collisionless bulk 
plasma, such as finite ion Larmor radii (FLR), finite parallel electric field (SEh), and 
kinetic ion compression (SG{). 

In the present paper, we exclude kinetic compression and construct a fluid limit, 
which we call "FLR MHD." This model is not self-consistent, but it is useful to isolate 
the effect of FLR terms on ideal MHD aTAEs. The questions dealt with in the present 
paper may be summarized as follows: 

(i) What are aTAEs? (a review with detailed discussion of bound state formation 
in ideal MHD toroidal plasma) 

(ii) How are aTAE affected by FLR corrections? (analysis of the Schrodingcr 
potential, mode structures and eigenvalues; in particular, frequency shifts) 

(iii) How is the propagating component of the mode structure affected by FLR 
corrections? (peculiar features of shear Alfven continuum waves in FLR MHD) 

The results of item (ii) will enable us to identify the branches of ion-temperature- 
gradient (ITG)-driven Alfvenic instabilities observed in gyrokinetic simulations solving 
the self-consistent model equations. This is done in the companion paper [? ]. 

The FLR MHD model and the numerical methods used are described in section 
2. In section 3, we review the physical origin and basic features of ideal MHD aTAEs. 
In section 4, we show how the structure of the effective Schrodingcr potential in 
the FLR MHD model varies with the mode frequency and demonstrate that, in 
cases of interest, the dominant bound state component of an aTAE cigenfunction 
remains essentially unchanged by FLR corrections. The eigcnfrcqucncies of two aTAE 
branches are inspected in section 5, where we analyze their dependence on a (pressure 
gradient), k$p c i (poloidal wavenumber x ion Larmor radius), and r]i (ratio of density 
and temperature gradient scale lengths). It is also shown that the magnitude of 
the diamagnetic frequency shift depends on the mode structure. In section 6, the 
results are summarized and conclusions are drawn. The Appendix contains a detailed 
discussion of how (propagating) continuum waves are modified by FLR effects. 

With our choice of terminology, we propose to use the term "toroidal Alfven 
eigenmode" (TAE) in a more general sense than is usually done. Our concept of TAE 
encompasses all ideal MHD cigenmodes formed inside the quasi-periodic structure of 
the effective Schrodingcr potential along a flux tube; which cannot be found in the 
cylindrical limit, only in toroidal geometry. This incorporates effects due to finite 
aspect ratio (e = a/Ro), shaping of flux surfaces (cllipticity, triangularity, etc.), and 
distortions caused by the pressure gradient (a) . The "classical" TAE is associated with 
the £-dependence of the magnetic field strength, B{0) w Bq/(1 + ecos9). When the 
a-induced potential barriers dominate (as in the present work), we speak of aTAEs. 
Naturally, there are parameter regimes where no clear distinction can be made between 
the "classical" TAE and aTAE: typically, when e and a have similar values. More 
generally, one can refer to Alfven cigenmodes (AE) as discrete bound states due to 
equilibrium non-uniformities causing poloidal symmetry breaking [? ] . 
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2. Physical model and numerical method 

The present work is part of an effort to study and understand the dynamics of drift 
Alfven ballooning modes in a wide range of frequencies and wavelengths and their 
kinetic excitation in a high-/? plasma with toroidal geometry. The full model is 
given by the one-dimensional linear gyrokinetic equations based on the derivation 
by Chen & Hasegawa [? ]. There, the electromagnetic field perturbations are 
described by gyrokinetic Maxwell equations in terms of the magnetic flux function 
Sip, the electrostatic potential 8<f), and the magnetic compression 8Bn . The linear 
gyrokinetic Vlasov equation governs the evolution of the fluctuating part, Sf s , of the 
total distribution function, f s = fo s + <5/ s , where the subscript s labels the particle 
species (not to be confused with the magnetic shear). 

Here, we utilize a reduced model which is obtained from [? ] (and more recent 
formulations in [? ? ]) by excluding energetic ions, ignoring kinetic thermal ion 
compression, including magnetic compression only at lowest order, and assuming a 
Maxwellian equilibrium distribution, /o s • We also ignore the variation of the magnetic 
field strength along a flux tube by letting B = Bo, so there are no magnetically 
trapped particles and no toroidicity-induced gap in the Alfven continuum. However, 
magnetic curvature and V-B drifts are properly accounted for. The reduced equations 
are presented in section 2.1, where the connection with the original gyrokinetic model 
is shown and the physical meaning of individual terms is given. In section 2.2, the 
three equations are combined into one single SAW equation, written in the form of a 
linear time-independent Schrodingcr equation. This equation is solved with a standard 
shooting method, the appropriate boundary condition for which is given in section 2.3. 



2.1. FLR MHD equations for shear Alfven waves 



With the use of the ballooning formalism [?????], the poloidal angle coordinate, 
$ G [0,27r], is mapped onto an infinite covering space, 9 G (—00,00). The coordinate 
9 effectively measures the distance along the field line normalized by qRo, with Rq 
being the major radius of the magnetic axis and q the safety factor measuring the 
average field line pitch. Thus, the connection length in 9 is 2tt. For convenience, we 
write the equations in Laplace-transformed form (dt — ► —iu, with 10 = lo t + i"f). We 
use SI units and define the following coefficients and parameters: 
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(vl + v?,)/2 is the kinetic energy and 
/i = v 2 \_/{2B) the magnetic moment. aj* s is the diamagnetic frequency associated 
with the density gradient, lo^s the magnetic drift frequency, and lo c0s the cyclotron 
frequency at the magnetic axis. The quantities /, g and h describe the flux 
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tube geometry in the shifted-circlc model equilibrium in terms of the parameters 
a = —q 2 R ft' and s = rq' /q, where (3 = 2{i P/B 2 is the ratio between thermal and 
magnetic pressure, and the prime denotes a radial derivative d/dr. Apart from a 
brief review of the effect of nonzero 9k (variable radial envelope) in section 3, we 
consider the case 9 k — 0. Here and in the following, the attributes parallel (||) and 
perpendicular (_L) refer to the direction relative to the equilibrium magnetic held, the 
toroidal component of which is assumed to be dominant as in typical tokamaks. 

Following standard procedures, the adiabatic and convective responses are 
separated through the substitution 

Sf. = ( H + ^J 5^e iL A + : (2) 



T 

S 1 .S 



where terms involving d^fos are omitted since the equilibrium distribution, /o s = 
n s0 (27r7^) -3 / 2 exp(— £ /T s ), is isotropic. J n (A s ) is the Bessel function of order n 
introduced by the gyroaverage. Lk = —k± ■ (b x v±)/uj c o s is the generator of 
the coordinate transformation between guiding center and particle variables, and 
e iLk turns into another Jo when gyroaveraged. The quantity SG S captures the 
compressional part of the non-adiabatic component of the particle response (in short, 
kinetic compression). Electrons are approximated as a massless fluid, so that 8G C = 0. 
The evolution of kinetic ion compression, SGi, is governed by the gyrokinetic equation 

[?]■ 

In the present study, we wish to ignore kinetic compression, so we let SGi = 0. 
This yields the following fluid equation, which we refer to as FLR MHD model: 
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where (...) = J d 3 v = J2a I d£ / d/iS/|«||| is the velocity space integral. Equation (3) 
is the so-called vorticity equation, which is obtained through the combination of the 
parallel Ampere's law with the continuity equation [? ? ]. Its individual terms 
describe field line bending (FLB), inertia (with FLR), MHD particle compression 
(MPC), and magnetic field compression (MFC). Both MPC and MFC arc static 
compression effects associated with toroidal curvature and finite (3. Equation (4) 
is the quasi-neutrality condition, which originally read J2 S ( e sSf s ) = 0. Equation (5) 
is the perpendicular Ampere's law reduced to an equation for the total (magnetic plus 
thermal) pressure balance for convective displacements. Note that oj6B\\ = —Cl p Sip 
effectively eliminates Q p (high-/? correction to the curvature drift) from Wdi at the 
order considered and reduces the drift-kinetic part of the "MPC+MFC" term to 
-2Q.pQ.Jip = agkl/(qR ) 2 8i! ("ideal MHD ballooning term"). 

Strictly speaking, the kinetic compression terms neglected here have a non-zero 
MHD limit if one goes to a regime where \ui\ ^> \u><n\ and \uj\ |fc||fi||. However, in 
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the parameter regime we are dealing with, these conditions are only satisfied for the 
low-energy part of the ion distribution. Thus, one cannot integrate over the entire 
velocity space as is usually done to compute the fluid limit (e.g., when estimating 
the width of the kinetic thermal ion gap [? ]); otherwise, the estimate will be wrong 
and one may as well neglect it altogether, as we choose to do here. The FLR MHD 
model may consequently be regarded as "incompressible FLR MHD," in the sense 
that the ion sound branch, compressibility associated with finite aspect ratio (trapped 
particles), and wave-particle interactions are not included. 



2. 2. FLR MHD equation in Schrodinger form 



In order to write the equations in dimensionlcss form, we employ the following 
normalizations : 
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where wao = v ao/(qRo) is the Alfven frequency. In the following, the hats will be 
neglected on all quantities except koi- 

The vorticity equation (3) is combined with equations (4) and (5), and brought 
into Schrbdinger-form through the substitution S^s = \/JSip, which yields 

- V cS (9)S^s = 0. (7) 

Here, J*" stands for d 2 (<5\I> s ) / dd 2 . The effective Schrodinger potential, 

V eS = V + V m , u + V m . T + K,flr, (8) 
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The functions To, Ai, Ti and T2 K are defined as (cf. equation (2.22) of [? ]) 
r = e -<%, A x = 1 + (h/Io - l)6i/2, 

Ti = 1 + {h/h - l)Vih, T 2K = A x + [1 + (3/i/Jo - 4)6i/2 - (h/I - l)bf] m; (13) 

with Ik{bi) being the modified Bcssel function fc-th order. 

The terms S^S" — V8^ s arc obtained by combining the "FLB" term with the drift- 
kinetic part of the "MPC+MFC" term (ideal MHD ballooning term) from equation 
(3). The ideal MHD potential V is determined by the equilibrium field geometry 
parametrized by s and a. The "inertia" term of equation (3) is now —(V mU) + V mtT ), 
where V myUJ consists of ideal MHD inertia, u> 2 , plus FLR corrections, and V myT captures 
the effect of non-zero parallel electric field, SEn = —dg(5<fi — Sip). The subscript "m" 
stands for mass (ion inertia). The FLR correction of the "MPC+MFC" terms is 
captured by V^flr- 
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2.3. Outgoing boundary condition 

At the outer boundary 9 = # max of the shooting range (in the asymmetric case, 9^ ^ 0, 
on both sides of the domain, 9 = ±# max ), we apply an outgoing boundary condition 
by matching the numerical solution to an analytical solution. The latter is obtained 
by solving equation (7) in the limit |fcoiS#| 3> 1, where it asymptotically approaches 

X(u) rT _,. . sin( 

(k is9) 2 k Qi s6 

The coefficients 



= S < + IT^L 5 ** + Y(")^M B , (\9\ » l/\k 0i s\) • (14) 



X(uS) = (oj + T ei LJ*i)(uj - w*i)/(l + r ci ), Y(w) = (u — w*pi)2qv t i (15) 

originate from — (V mi u+Vm iT ) and — V k .flr, respectively As is shown in the Appendix, 
the solution of equation (14) which satisfies the causality constraint for an outgoing 
group velocity, v g /9 > 0, is 
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where a = Re{C}/|Re{C}|, C = y/AX + 2Y 2 - 1 and ^ is a constant. Formally, we 
write 8\V' s /5ty s = ikii. By matching the analytic solution given by equation (16) to the 
numerical solution for equation (7), one ensures that the waves which travel through 
the entire simulation domain, are not reflected at the artificial boundaries. 

Other than this, equation (14) and its solution has no physical meaning, because 
it only describes the balance between magnetic induction and the ion polarization 
current. The latter vanishes in the limit |fco;s#| ^ 1; where the radial width of the 
perturbation becomes much smaller than the ion Larmor radius. Physically, magnetic 
induction will then be balanced by electron inertia or the wave is dissipated by 
collisions. These effects are not included in equation (7) and the underlying gyrokinctic 
Maxwell- Vlasov equations, and they are not important for the study of oTAEs which 
are localized near \9\ ~ 0{1). oTAEs may lose energy to the outgoing waves, which 
represents the physical effect of continuum damping, but are not otherwise affected 
by their presence. The physics of large \9\ become important only near marginal ideal 
MHD ballooning stability, where the mode structure becomes broad. 



2-4- Parameter values 

In the present paper, we analyze two cases: one with lower magnetic shear (s = 0.4) 
and one with higher magnetic shear (s = 1.0). The default parameters are listed in 
table 1. They were adopted (with minor changes) from two earlier linear gyrokinctic 
studies of Alfvenic ITG instabilities [? ? ] . 

Note that is is not meaningful to simply compare the values of parameters (s, a) of 
the shifted-circlc equilibrium model used here with parameter values in experimental 
configurations with non-circular flux surfaces. The location of the stability boundaries 
of the ideally unstable domain and the distribution of the oTAE bands in the s-a plane 
are sensitive to the flux surface geometry, so it is likely that a given region in parameter 
space of the simple model may correspond to another region of the parameter space 
in a more realistic geometry. A meaningful comparison would need to consider the 
shape of the Schrodinger potential at different points in parameter space, rather than 
the values of parameters like s and a. This may be the subject of a future study. 
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Table 1. Default physical parameters in the two cases considered in this paper: 
one with lower magnetic shear (s = 0.4) and one with higher shear (s = 1.0). 
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ss a/80 


Dong et at. [? ] 



3. Review of ideal MHD Alfven eigenmodes 

Many qualitative properties of ideal MHD aTAEs are preserved in FLR MHD. Thus, 
it is instructive to review the physics underlying aTAEs in the simpler ideal MHD 
limit, where the SAW Schrodinger equation (7) reduces to 

= J< - V5^ s + u 2 S^ s = [(fS^'Y + agSiP + u 2 f5ijj] jyff. (17) 

In section 3.1, we elucidate the origin of the a-induced potential barriers. The types of 
discrete shear Alfven eigenmodes are introduced in section 3.2: purely growing ideal 
MHD ballooning modes and damped oscillatory solutions, now called aTAEs. The 
latter are described in detail in section 3.3. Section 3.4 contains remarks regarding 
the effect of the radial envelope and global trapping. 

3.1. Equilibrium magnetic field geometry 

We employ the s-a model equilibrium for simplicity and comparability with earlier 
works (the results of which we explain in the companion paper [? ]). The model is not 
rigorously valid for a > 1, but even outside its regime of validity it captures essential 
qualitative properties of more accurate models [? ? ]. The most prominent feature 
predicted by the s-a model is the second MHD ballooning stable domain. 

The model is useful because it captures two essential ingredients: toroidicity and 
magnetic shear. The resulting inhomogeneities and asymmetries give rise to discrete 
shear Alfven eigenmodes which do not exist in the cylindrical limit. Writing the 
vorticity equation in Schrodinger form, such as equation (17), offers an illuminating 
way to express these physics in terms of familiar quantum mechanics concepts: 
toroidicity and magnetic shear introduce quasi-periodic barriers and wells in the 
effective Schrodinger potential along a flux tube. As a consequence, the spectrum 
of ideal MHD shear Alfven waves in toroidal geometry consists of continua (waves 
propagating along a flux tube) and discrete eigenmodes (standing waves satisfying 
quantization conditions). Further geometric complications introduce additional gaps 
in the continuum and corresponding branches of discrete modes at various frequencies. 

Speaking in simple terms, in the s-a model, the periodicity of the barriers is 
controlled by the magnetic shear s and their height by the pressure gradient parameter 
a. For the discrete eigenmodes studied here, the dominant effect of toroidicity is the 
one that enters via a. We look at a regime where the a-induced potential barriers are 
significantly larger than those induced by the inverse aspect ratio a/ Rq, so we neglect 
the variation of B and R along a flux tube (along with all the effects of non-circular 
flux surfaces) and let B{6) = B and R(6) = R (constant). 

The potential structure may be traced back to its physical origin as follows. The 
function g = cos 8 + h sin 8 appearing in the ideal MHD ballooning term represents 
the sum of normal and geodesic curvatures. The factor h reflects the fact that the 
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modulation of geodesic curvature is intrinsically coupled to the local magnetic shear, 
which is measured by the quantity h' = s — a cos 9. Noting that the term —a cos 9 in h' 
has the form of the Pfirsch-Schliiter current [? ] , it is clear that the s-a model describes 
the modulation of local magnetic shear and geodesic curvature which is produced by 
currents that balance the charge separation induced by toroidal curvature drift. 

Taking into account these distortions of the flux tube geometry, the Schrodingcr 
potential V(#|s, a) measures how much energy is needed to bend a field line at a given 
location 9. From the functional form of V defined in equation (9), the numerator 
of which is (s — a cos 9) 2 — fa cos 9 (with / > 1), it can be seen that potential wells 
(V < 0) are found at those locations where the average magnetic shear s is canceled by 
the distortions caused by the Pfirsch-Schliiter current. Potential barriers (V large and 
positive) are found in those regions where field line bending is most strongly inhibited 
due to increased local shear. 

3.2. Branches of ideal MHD Alfven eigenmodes 

The potential V measures how much energy is needed to bend magnetic field lines. 
Increasing field line bending energy makes V larger, while increasing free expansion 
energy associated with the pressure gradient a makes V smaller. In the potential wells 
where V is negative, free expansion energy exceeds field line bending energy locally. 
This occurs in regions with unfavorable curvature with respect to the direction of the 
pressure gradient when the local magnetic shear, \h'\, is sufficiently small. 

A shear Alfven wave travelling along a field line and encountering a barrier will 
be reflected, at least partially. Thus, wave energy tends to be trapped between 
the barriers and standing waves can form whenever the wave function satisfies the 
quantization condition imposed by the barriers. There are two types of discrete 
solutions: purely growing eigenmodes and damped oscillatory eigenmodes. The 
damping of the latter is due to the finite width and finite height of the potential 
barriers, so there is either tunneling or direct coupling to waves propagating outward 
(corresponding to the continuous spectrum). The purely growing solutions are those 
which can tap a sufficient amount of free expansion energy to overcompensate field 
line bending energy The Schrodingcr potential V consists of an infinite number of 
potential barriers and wells, so the s-a plane contains overlapping bands of discrete 
shear Alfven eigenmodes localized in one or several potential wells. In the following, 
we shall consider some concrete examples. 

Purely growing solutions are known as ideal MHD ballooning modes (BM). Two 
bands of purely growing ballooning modes, having different radial envelopes 9k = 
and 9k = n, are shown in figure 1(a). Unstable bands with intermediate values 
< \9k\ < 7r are located in the region spanned by the two bands shown [? ]. In 
figure 1(b), the growth rates in these two bands are plotted as functions of a for 
s = 1.0. The corresponding mode structures and potentials are plotted in figure l(i) 
and (ii). The potential structure is symmetric around 9 = 9k for integer values of 
9k /it and asymmetric otherwise. 

For a given value of s, the ballooning branch appears at a cr it.i(s) [~ 0.6 for 
band (i) in figure 1(b)] and merges back into the accumulation point u> — at the 
second stability boundary, a cl it,2(s) [~ 2.6 for band (i) in figure 1(b)], as field line 
bending again balances free expansion energy. As a increases past this point, the 
purely growing ballooning branch is replaced by an oscillatory bound state localized 
in the same potential well, — it < 9 < tt. This oTAE branch is labeled (iii) in figure 




Figure 1. Ideal MHD shear Alfven cigenmodes in the s-a plane, (a): Examples 
of ideal MHD ballooning unstable bands (BM=ballooning mode), for 9^=0 and 
@k = n - (b) : a-dependence of the growth rates 7 of the BM bands 9^=0 (i) and 
&k = n (ii) f° r s = 1-0' The vertical dotted lines indicate the first and second 
stability boundaries for 9 k = 0, which are labeled a cr it 1 and a cr it 2> respectively, 
(c): ct-dependence of the real frequencies aj r of aTAE(j,p) ground states (p = 0) 
for 9^ = 0. ctTAEs localized in the central (iii) and in the second potential 
well (iv) are shown. a num is the numerical threshold above which band (iv) 
could be identified. On the right-hand side, diagrams (i)-(iv) show the typical 
structure of the Schrodinger potential V(&) and cigenmodes |<5\]> s (0)| in covering 
space —00 < 9 < 00 for the bands (i)-(iv) shown in (b) and (c). 



1(c), where the eigenfrequency is plotted as a function of a. As a increases further, 
another aTAE localized in the second potential well, 7r < \9\ < 37r, is found. In figure 
1(c), this branch is shown for a > 7 and labeled (iv). The typical mode structures are 
plotted in figure 1 (iii) and (iv) along with the corresponding potential V. 

The range of a values scanned exceeds by far the range of validity of the s-a 
model and probably goes far beyond what is realistically achievable in a tokamak- 
type plasma. However, this exaggeration is necessary to reveal the band structure of 
the s-a plane and is useful to study the properties of the shear Alfven eigenmodes 
which populate it. 

3.3. Characterization of a-induced toroidal Alfven eigenmodes (aTAE) 

In general, aTAEs consist of a propagating component, which belongs to the Alfven 
continuum, and a standing wave component trapped between potential barriers 
induced by the pressure gradient measured by a. aTAEs can exist as bound or 
unbound states. Mathematically, the distinction is made on the basis of the location 
of turning points in the complex 6 = 6 T +i6{ plane of the Stokes diagram [? ] (cf. figure 
9 in [? ]). For a bound state, the turning points associated with the trapping potential 
barriers lie on the 9 T axis, whereas turning points of unbound states have 8i 7^ 0. In 
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other words, for unbound states, the reflected component is coupled to the continuum 
directly (and, thus, suffers strong damping), whereas bound states transmit energy to 
the continuum only via barrier tunneling. 

Since there is an infinite number of potential barriers, there exists a dense 
spectrum of aTAEs. However, in practice, only a few modes are distinguishable from 
the continuum, and these are found only for sufficiently large values of a; i.e., primarily 
in the second MHD stable domain. Following the notation introduced in [? ] , aTAE 
branches are denoted as aTAE(j, p) . The label j > 1 identifies the potential well where 
the mode peaks, with potential barriers approximately located at 9 ~ ±7r(2|j — 1| — 1) 
and 9 ~ ±7r(2j — 1). The label p > counts the number of zeros the mode structure 
has in potential well j and corresponds to the energy quantum number, with p = 
being the ground state. The eigenfrequency is denoted by w^j^ E . The two branches 
in figure 1(c) are oTAE(l,0) [branch (iii)] and aTAE(2,0) [branch (iv)]. 

Generally speaking, reflections off a-induced potential barriers occur for any a > 
0; first, they form unbound states when a is small; later, when a is sufficiently large, 
these turn into bound states. The aTAE(l, 0) branch is a special case: it is replaced 
by the unstable ideal MHD ballooning branch in the region a cr it,i < a < a cr ;t.2- For 
sufficiently large values of s > 0.5, the oTAE(l,0) branch reappears at a = a cr it.2 in 
the form of a strongly bound state. 

One may imagine similar bands of ideal MHD ballooning modes in potential wells 
with j > 1, but no such modes were found in the s-a model equilibrium. Thus, as 
a is increased, aTAE branches with j > 1 (as well as states with j = 1 and p > 0) 
smoothly transform from an unbound to a weakly bound and, eventually, strongly 
bound state. When solving equation (7) with the shooting method, a has to exceed 
a certain numerical threshold, a num , beyond which an aTAE is sufficiently deeply 
trapped to be clearly identifiable as a distinct discrete eigenmode. In the case of the 
oTAE(2,0) branch in figure 1(c), we have a num ~ 7. 

As aTAEs become strongly bound ( "quasi-marginally stable") with increasing a, 
they attain high frequencies comparable to the Alfven frequency. The two examples 
(iii) and (iv) in figure 1(c) illustrate this fact. Naturally, modes with larger energy 
quantum number p tend to reach higher frequencies due to stronger field line bending. 
In the s-a plane, branches with different j and p can become degenerate by crossing 
or merging (not shown). For instance, the even (1,0) branch merges with the odd 
(1,1) branch when a becomes sufficiently large, because a potential barrier appears 
at 9 = 0. 

3.4- Effect of the radial envelope 

Our use of the term "bound state" only refers to the trapping of waves along a flux 
tube. An inspection of the dependence of the eigenfrequencies on the minor radial 
coordinate, r, and the wave number of the radial envelope, 9k, is necessary to determine 
whether a given mode is a globally bound state or whether it propagates radially [? 
]. A necessary (not sufficient) condition for the existence of a globally bound state 
is duj r /d9k = 0. This is only the case when the potential V is symmetric around 
the location where the mode is locally trapped. For aTAEs trapped in the central 
potential well (j = 1) this is the case for 0^ = 0, and for the second well (j = 2) the 
condition is satisfied for 9k = n. aTAEs with j > 2 cannot form globally bound states 
because du> r /d9k ^ for any 9k (note that 9^ — * 9k ± 27r implies \j — 1| — ► \j — l^f 1|). 
A higher excitation threshold due to radial propagation may be expected for these 
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modes. Further global analysis goes beyond the scope of the present paper. In the 
remainder of this work, we let Ok = (flat radial envelope) in order to be able to make 
comparisons with earlier local gyrokinetic studies. 

4. Effect of FLR on the Schrodinger potential and aTAE eigenfunctions 

In this section, we examine how FLR effects modify the effective Schrodinger potential 
(section 4.1) and aTAE mode structures (section 4.2). It is found that the structure 
of V e s(9) is similar to the ideal MHD potential V(6) when uj w oj* P i; otherwise, the 
potential structure is modified by the term Vflr,k given in equation (12). In the 
eigenmode structure, no significant change is observed in the shape of the dominant 
bound state component. However, FLR terms may reduce the damping rate by moving 
additional turning points to the 8 r axis of the Stokes diagram and strongly modify the 
waveform of the outward propagating component. 

4-1. Shape of the Schrodinger potential 

The shape of the ideal MHD potential V(6\s,a) is uniquely determined by the 
equilibrium geometry parametrized by s and a. When the FLR terms are included in 
equation (7), the effective potential acquires dependencies on many other parameters: 
Vos (0\u), s, a, rji, h x , e n , q, Vt\) [cf. equation (8)]. In particular, note that the dependence 
on uj and fcoi implies that modes with different frequencies and wavelengths observe 
different potential structures. This is illustrated in figure 2, where the effective 
potential V g is shown for several values of u> and compared with the ideal MHD 
potential V. Despite this complexity, we can make several general remarks: 




Figure 3. Mode structure of a quasi-marginally stable »TAE(1,0) for the case 
s = 1.0 (cf. table 1) with a = 3.9. Results are shown for (a) the ideal MHD 
limit [equation (17)] and (b) for FLR MHD [equation (7)], where uj tpi = 0.88. 
The respective effective Schrodingcr potential V c ff ( — ■ — ) is plotted along with 

the real ( ) and the imaginary component of <5*l/ s ( ). Panel (c) shows 

the Stokes diagram with local anti-Stokes lines ( ) for the ideal MHD (top 

half) and FLR MHD case (bottom half). In addition, global anti-Stokes lines 

which lie on the 8 T axis (bold ) and the associated turning points (• ) are 

indicated. Anti-Stokes lines represent paths in the complex 8 = 8 r + i8[ plane 
along which J dflQ 1 / 2 is imaginary (here, Q = —V e g, u = ui t ); i.e., where the 
cikonal approximation yields purely oscillatory solutions [? ]. 



• FLR terms may modify the shape of the bottom of a potential well (modifying the 
ballooning stability limits) and the shape of minor potential barriers (modifying 
their wave trapping properties). In the absence of kinetic compression, the effect 
on ballooning modes is found to be a stabilizing one when r/j > [? ]. 

• The potential V e g is closest to its ideal-MHD limit V when to ~ as can be 
seen in figure 2(c). This is particularly true for \0\ ^> l/\ko[s\, as can be seen from 
the form of Y{ui) in equation (15) [which orginates from Vflr.k in equation (12). 

• FLR terms add offsets to the cigcnfrcqucncy, such as the well-known diamagnetic 
frequency shift. As is shown in section 5.3 below, the effective value of this offset 
can vary from one potential well to the next, so the frequency shifts depend on 
the mode structure (i.e., the location of the dominant bound state component). 
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4-2. aTAE eigenf unctions 

We analyze the effect of FLR terms on the mode structure of aTAEs using the 
examples in figures 3 and 4, where an aTAE(l,0) for higher shear (s = 1.0) and 
an aTAE(2,0) for lower shear (s = 0.4) are shown, respectively. The figures show (a) 
the ideal MHD result, (b) the FLR MHD result, and (c) the Stokes diagrams for the 
solutions in (a) and (b). As the Stokes diagrams show, the eigenmodes in both cases 
are strictly bound states; in fact, the aTAE(l,0) in figure 3 is already in the strongly 
bound (quasi-marginally stable) regime. The effective potential, V c g, and eigenvalues, 
u> = ui T + ij, are also shown in panels (a) and (b) of both figures 3 and 4. 

No significant change in the structure of the dominant bound state components 
[identified by the labels (j,p)] is observed. Hence, we may assume that the increase in 
the eigenfrequency w r in figures 3 and 4 is not due to a change in field line bending, 
but only an FLR-induccd offset, which is discussed in section 5 below. 

Energy which tunnels through the potential barriers is carried away by the 
outgoing continuum wave which matches the bound state's frequency (resonant 
absorption). Most of the difference in the mode structures shown in figures 3 and 
4 lies in the waveform of this propagating component: 

• Panel (a): Ideal MHD continuum waves are harmonic wave solutions obtained in 
the large- |s#| limit of the ideal MHD equation (17). The finite damping rate of an 
aTAE is reflected in an exponentially growing amplitude of the form cxp(— C7|0|), 
where C is a real positive number. 

• Panel (b): The FLR continuum waves are governed by equation (14) obtained 
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for \9\ l/|fcoi s l °f the FLR MHD model (7). Distinctive features arc an 
algebraically increasing amplitude (roughly oc l^l 1 / 2 ) and an increasing phase 
velocity and wavelength (both oc \6\) corresponding to the decay of the ion 
polarization current. A more detailed discussion is given in the Appendix. 

These differences have to be taken into account when imposing boundary conditions 
in numerical codes used to study weakly bound aTAEs such as those in figure 4 (more 
about this in the Appendix). This is necessary to accurately compute the continuum 
damping, which may contribute to the instability threshold in the presence of wave- 
particle interactions, as is discussed in the companion paper [? ]. The eigenfrcqucncics 
are almost insensitive to the treatment of the outgoing continuum wave, because they 
are determined between the turning points associated with the primary potential well, 
j, where the main component of the aTAE(j,p) is trapped. 

The bold solid lines in panels (a) and (b) of figures 3 and 4 represent the envelope 
of the mode structure. The modulation of this envelope indicates additional trapping 
of wave energy in secondary potential wells with j + j + 2, etc. In the case with higher 
shear, shown in figure 3, we see that FLR terms significantly enhance the amount of 
trapping in secondary wells. The Stokes diagram in figure 3(c) reveals that, in the 
FLR MHD case, additional turning points are present on the 6* r axis. This is consistent 
with the fact that the damping rate drops by one order of magnitude: from —7 w 10 -4 
in ideal MHD [figure 3(a)] to -7 « 1CT 5 in FLR MHD [figure 3(b)]. 

In contrast, in the case with lower shear, shown in figure 4, there is no significant 
change in the amount of trapping in secondary wells; the envelope \*5f s \, the Stokes 
diagram, and the damping rates are similar in ideal MHD (a) and FLR MHD (b). 
According to the Stokes diagram in figure 4(c), both the ideal MHD and the FLR 
MHD eigenmode may be viewed as a combination of a primary (2,0) and a secondary 
(3,0) bound state component. 

5. Effect of FLR on aTAE eigenfrequencies 

In this section, we inspect how FLR effects modify the eigenfrequencies of aTAEs. 
Eventually, the eigenfrequencies obtained here with the FLR MHD model (7) will 
allow us to identify Alfvenic ITG instabilities seen in gyrokinetic simulation of the 
second ballooning stable domain [? ] as aTAEs, which is done in a companion paper 
[? ]. With this motivation in mind, we study parameter scans with respect to a, 
rji and fcoi in sections 5.1 and 5.2. In section 5.3 it is shown that, when 77; > 0, the 
magnitude of the diamagnetic frequency shift depends on the mode structure; i.e., the 
location of the dominant bound state component. 

5.1. a- dependence of the eigenfrequency 

Figure 5 shows the frequencies of two branches, aTAE(l, 0) ( — • — ) and aTAE(2, 0) 
( -), as functions of the normalized pressure gradient a. Ideal MHD results are 
shown in figure 5(a) and (b) for s = 1.0 and s = 0.4, respectively. As discussed in 
section 3, our method does not allow to follow the aTAE(2,0) branch below a certain 
numerical threshold, a 11U m- Typical mode structures | <5\l/ s | found on the two branches 
are shown in panels (i) and (ii). 

Corresponding results obtained with the FLR MHD model (7) are shown in figure 
5(c) and (d) where, in addition to the aTAE branches, the frequency of the FLR- 
modificd ballooning mode (FLR BM) is plotted as well (• ). For comparison, the 
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Figure 5. a-dcpendence of oTAE eigenfrequencies in (a,b) ideal MHD and (c,d) 
FLR MHD for two values of magnetic shear: (a,c) s = 0.4 and (b,d) s = 1.0 

(cf. table 1). In addition to the aTAE branches (1,0) ( — ■ — ) and (2,0) ( ), 

the frequency of the ideally unstable FLR ballooning mode (FLR BM) is shown in 
(c) and (d) (• ). For comparison, the diamagnetic frequency ui tp i is plotted in (c) 

and (d) ( ). The ballooning stability boundaries ( ) are labeled a cr it,i 

and a cr it 2i the numerical threshold for identifying an aTAE branch is labeled 
ci llnm) and the label ao indicates the point at which the curve uj T (a) exhibits a 
kink. Mode structure examples are shown in panels (i)-(iv). 



diamagnetic frequency uj^a is also shown ( ). Typical mode structures \8^ s \ found 

on the BM and aTAE(2,0) branch are shown in panels (iii) and (iv). Note that the 
numerical thresholds a num (s) are reduced in the FLR MHD case, so we are able to 
identify distinct aTAE(2,0) bound states at lower values of a than in the ideal MHD 
limit. All aTAE solutions shown in figure 5 are strictly bound states; with turning 
points on the 8 T axis of the Stokes diagram (not shown). 

When comparing the FLR MHD case [panels (c) and (d)] with the ideal MHD 
limit [panels (a) and (b)] in figure 5, several interesting features can be observed: 

• As expected, the ballooning unstable domain shrinks and the eigenfrequencies 
undergo an up-shift. The frequency shift is primarily due to the diamagnetic 
frequency o;* p i, with a significant contribution from Vflr.t due to finite 5E\\ 
(finite tJ{). When trying to quantify the frequency shift, we find that its value 
depends on the mode structure. This phenomenon is discussed in section 5.3 
below. 

• Near the second FLR MHD ballooning stability boundary, the frequency of 
the FLR aTAE(l,0) branch drops well below cj* p i, which is known to be the 
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Figure 6. FLR aTAE(2,0) cigcnfrcqucncies, lj t (O ), and damping rates, —7 
(x), as functions of a for (a) s = 1.0 and (b) s = 0.4 [cf. figure 5(c) and (d)]. The 
location where the the slope of the curve u! T (a) changes is labeled oq ( ). 
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Figure 7. Eigenvalues and cigenfunctions of an «TAE(2,0) as functions of rji (left 
panels) and fcgi (right panels) for the case with s = 1.0 (cf. table 1). The (negative) 
growth rates are shown in (a,c), the frequencies in (b,d), and mode structures in 

(e,f). For comparison, the diamagnctic frequency is plotted in (b) and (c) ( ). 

The shaded areas in (e) and (f) highlight the bound state component. 



continuum accumulation point of the diamagnctic gap in lowest-order FLR 
theories [? ]. This observation is also explained in section 5.3 below. 

• In the FLR MHD case, the curves w r (a) for the aTAE(2, 0) branch exhibit a kink 
at a location labeled a — uq. For a < ao, the frequency approximately scales like 
uj r cx a 1 / 2 ] whereas a steep rise in the frequency occurs for a > ao- A comparison 
with the behavior of the damping rate —7(a) in figure 6 shows that, at least in 
the case with higher shear, ao coincides with the threshold beyond which the 
bound state may be regarded as deeply trapped ( "quasi-marginally stable" ) . 

5.2. 7]i- and koi- dependence of the eigenfrequency 

Using the data point a = 6.7 on the FLR o;TAE(2,0) branch in figure 5(c) as a starting 
point, we vary the parameters rji and koi. The results are shown in figure 7. 

In the entire parameter range scanned in figure 7(a)-(d), the modes are strictly 
bound states and are subject to damping due to barrier tunneling; with —j/uj r ~ 10 -2 . 
A sudden increase in the damping rate, such as near koi ~ 0.4 in figure 7(c), typically 
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Figure 8. Modc-structure-dependonce of the diamagnctic frequency shift. The 
examples are based on the case s = 1.0 in table 1 in the limit = 0. In panels 
(a)-(d), results are shown for an oTAE(l,0) at a = 2.9. In panels (e)-(h), results 
are shown for an aTAE(2,0) at a = 8.9. Panels (a) and (e) show the spatial 
dependence of the coefficient D u = (1 — ro)/&oi, and (b) and (f) that of the 
coefficient D* = (1 — r Ti)/[6 i(l + rfr)], for fc oi = 0.01 (— ■ — ) and k oi = 0.212 

( ). For fcoi = 0.212 [where uj tp i = 0A5y/a\, (c) and (g) show the ratio 

mi* = D*/D^ ( ), its asymptotic limit for \9\ — » oo. which equals 1/(1 + Vi) 

( ), and the normalized effective diamagnctic frequency w, a g = uj tp i/uj tp i 

( ) as seen by the respective eigenfunction. The eigenfunctions 8^/ s ( ) 

are plotted in (d) and (h) together with the corresponding ideal MHD potential 
V ( ). 

indicates a degeneracy with or a switching to an aTAE branch in a neighboring 
potential well, which is not necessarily evident in the eigenfrequency traces. In the 
case shown in figure 7(c), an analysis of the Stokes diagram (not shown) reveals that 
the dip in 7(fco0 near fcoi ~ 0-4 is correlated with the merging of a turning point with 
the 9 r axis near \9\ ~ 5ir; i.e., the formation of an additional bound state component 
(j,p) = (S,0). 

In figure 7(b) and (d), the frequencies are found to be near in most of the 
parameter range scanned. Although uj? ~ uj* p { has special implications for the shape of 
the effective potential [cf. section 4.1], we have not yet found any conclusive evidence 
showing that this is indeed a preferred eigenfrequency for aTAE in the regime a < olq . 

Note that the range of values of 771 and fcoi examined in figure 7 is relevant for 
Alfvenic ITG instabilities. Thus, the low damping rates and frequencies near the ion 
diamagnctic frequency suggest that these modes may be easily excited in the presence 
of ITG and kinetic compression, provided that uj r is also comparable to characteristic 
frequencies of particle motion. This is confirmed by gyrokinetic simulations [? ? ] . 

5.3. Mode-structure-dependence of the diamagnetic frequency shift 

Low-order FLR models based on a small-&i expansion, such as equation (A. 6) in the 
Appendix, predict that the diamagnctic frequency causes an up-shift of the entire 
SAW spectrum and a gap in the continuous spectrum. The frequency shift depends 



a TAE: I. Ideal MHD and FLR effects 



18 



only on tu« p i- When FLR effects arc fully retained, wc find that the frequency shifts 
also depend on the mode structure, because the terms responsible for the shifts are 
functions of 9, with a tendency to decay with increasing \9\. Typically, for modes 
localized in potential wells with j > 1, the full FLR model yields cigcnfrcquencies 
which are significantly lower than those obtained from low-order FLR models. Since 
the wave-particle resonance condition depends on the frequency, this effect is relevant 
for the resonant excitation of aTAEs [? ] . In this section, mode-structure-dependence 
of the diamagnetic frequency shift in the FLR MHD model (7) is demonstrated through 
comparison between an oTAE(l,0) and an oTAE(2,0). 

For simplicity, let us consider the case with SE^ = 0. In this limit, which is 
realized by letting = 0, the inertia term in equation (3) is 

-V m , u 6V a = (D^uj 2 - D*uuj t pi) (5* s . (18) 

After multiplication by the complex conjugate S 1 ^* and integration over 9, wc can 
define an effective diamagnetic frequency uupl as 

w{lo - ZU~) = u> (oj - uj*pi (D») 9 j (D u )y) ; (19) 
where the coefficients and D* are [cf. equation (10)] 

D u = (l-T )/h, D* = (l-r T 1 )/[(l + 771)61], (20) 
and the brackets represent a weighted average, 

= (j\eD^\6^ s A j {J\e\m s \ % if^ . (21) 

The limits of the integration interval, [0i,02]i correspond to the turning points. 

Two examples are presented in figure 8. Based on the case with s = 1.0 in table 
1, results are shown for a = 2.9 (a)-(d) and a = 8.9 (e)-(h), where deeply trapped 
aTAE(l,0) and oTAE(2,0) are considered, respectively. 

In figure 8(a) and (b), D u and -D* are plotted for fcoi = 0.01 and fcoi = 0.212, 
which clearly shows that FLR effects become important already at \9\/tt ~ 0(1). For 
fcoi = 0.212, the ratio = D^/D^ is plotted in figure 8(c) along with the normalized 
effective diamagnetic frequency, w t . c s = U^pl/u^pi, as seen by the aTAE(l,0) shown 
in figure 8(d). The results for the aTAE(2,0) in figure 8(e)-(h) are arranged in the 
same manner. 

It is found that, for the oTAE(l,0) in figure 8(d), the effective diamagnetic 
frequency amounts to 95% of ct>* p i = 0.76, so there is no significant change. In contrast, 
for the oTAE(2,0) in figure 8(h), the effective diamagnetic frequency is only 63% of 

= 1-33, so the associated frequency up-shift is correspondingly smaller. 

The effect is also visible in figures 5(c) and (d): at the second ballooning stability 
boundary, the transition between the FLR BM branch and the FLR aTAE(l,0) branch 
occurs at a frequency well below w* p i, close to w*; = w*j,i/(l + r/i). This may be 
readily understood by noting that, near marginal stability, the solutions acquire a two- 
scale structure, consisting of a localized peak [\ku\ ~ l/(qRo)] and a long-wavelength 
component (\ku \ — > 0) [? ]. The long- wavelength component dominates in the spatial 
average given by (21), so that a>* P i/a;*pi — > 1/(1 + 771 ) , which is the lower bound for 
the effective diamagnetic frequency. Obviously, this effect is only present for 771 > 0, 
since Ti(fii = 0) = 1 in equation (20). 



a TAE: I. Ideal MHD and FLR effects 



19 



6. Conclusion 

The properties of a-induced toroidal Alfven cigenmodes (aTAE) have been 
investigated in detail using a fluid model including full finitc-Larmor-radius (FLR) 
effects for thermal ions. The resulting changes in the effective Schrodinger potential, 
the eigenmodes and eigenvalues were described and explained. By isolating the physics 
of the fluid limit, the results reported in the present paper provide the foundation 
for understanding how non-resonant and resonant effects of kinetic ion compression 
further modify the eigenvalues and eigcnfunctions of aTAEs. 

aTAEs become quasi-marginally stable and acquire high frequencies to 3> Lu* p i 
when the normalized pressure gradient exceeds a branch-specific threshold, a > 
ao(j,p). In this regime, one may expect interactions with a population of energetic 
particles produced by auxiliary heating or nuclear fusion. Below this threshold, 
Q < &o(j,p), continuum damping is stronger but the eigenfrequency tends to be 
near the diamagnetic frequency lu^ and scales like a 1 / 2 . In this regime, aTAEs may 
be expected to be easily excited via interactions with thermal ions in the presence of 
an ion temperature gradient (ITG). Indeed, as is shown in detail in the companion 
paper [? ] , FLR- modified aTAEs are the normal modes which constitute the Alfvenic 
ITG instabilities first reported by Hirose et al. [? ] in the second ballooning stable 
domain. 
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Appendix: Outward propagating wave in FLR MHD 

This appendix summarizes peculiarities of the FLR MHD model which were 
encountered during our study. This model-specific information is included as a 
reference for researchers trying to reproduce our results and understand them in detail. 

Solution for the FLR continuum 

Equation (14) may be solved approximately with the ansatz 5^/ s =ip-\-ip sm.9, where 
ip is a long- wavelength component associated with the X/9 2 term and -0 is a short- 
wavelength component associated with the Y smQ/Q term. We assume the following 
formal ordering with respect to the small parameter 5 ~ l/|fcoiS#| <C 1: 



where the prime denotes a derivative with respect to 6. Separating the equations with 
respect to the long and short wavelengths, we obtain, at lowest order, 



i?/^~0(5), i>'/$~0{5), $/i)~0{5), \k oi s6\ ~ GOT 1 ); 



(A.l) 
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Figure Al. (a) Dispersion relation and (b) modulation of the algebraically 
growing amplitude for outgoing waves at large \0\. Results are shown for values of 
9 where the short-wavelength contribution <x Y cos 8 in equation (A. 3) vanishes 

( ), is maximal ( — • — ) and minimal ( ). The parameters for the case 

with s = 1.0 in table 1 are used, and we set 7 = —0.1, a = 6.3. For comparison, 
horizontal dotted lines in (a) and (b) indicate u)»p; and w f [, and the vertical dotted 
line in (b) marks the average algebraic growth associated with Im{fc|| } for the case 
7 = 0, which is l^ 1 / 2 . 



The solution for tjj has the form \k is9\( 1±iC ^ 2 with C 
so that 



[(4X + 2F 2 )/(fc 0iS ) 2 -l] 1/2 , 



1 1. 4-i — 



1 



Fsinl 
kois9 



11 26> 



±C- 



2Y cos 9 
k oi s 



(A.3) 



where kn = — i6^' s /S^> s and, typically, C 2 > 0. Equation (A.3) describes the 
propagating component of an eigenmodc which couples to the local shear Alfvcn 
continuum in FLR MHD, either through direct coupling or barrier tunneling. 

In order to determine the physically relevant sign in the exponent, we compute 
the group velocity u g = dw/dfcii. For simplicity, let us first consider a special type of 
solutions characterized by lo ps uj* p i. In this case, Y ~ and we have 



Re 



dfc|, 
duj r 



1 

±26 



— Re{C} 
dw r 



(A.4) 



so that 



v g P3 ±2|fcois|6»L>; (A.5) 

where D > is a constant. The outgoing wave has v g > for 9 > (and vice versa), 
so that the physically relevant solution to be matched at the boundary is the one with 
the positive sign in the exponent as in equation (16). 



The WKB dispersion relation 



k \\ = 



"II 



for any uj r is plotted in figure 



Al(a). The growth rate 7 enters mainly via the short- wavelength correction oc 
21 r cos6 , /(|fcois|6 l ) in equation (A.3), the real part of which vanishes when lo t = 
Although its magnitude is usually small, this short-wavelength correction can still be 
important for matching at the boundary because it determines the local slope of the 
wave function (unless u> T ~ 3> I7I). The short-wavelength modulation of the 
algebraic growth associated with Im{fc|i} is shown in figure Al(b). For 7 = 0, the 
average algebraic growth goes as l^ 1 / 2 , as can also be seen from equation (A.3). 
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Properties of the FLR continuum wave 

The properties of a full FLR model at large \9\ were previously analyzed by Connor 
et al. [? ] in a simpler geometry (sheared slab), which exhibits similar qualitative 
features as the toroidal flux tube model used here. In particular, since the inertia 
term in equation (14) disappears in the limit |fcoiS#| 3> 1, the full FLR model is very 
different from widely used lowest-order FLR models which are based on the assumption 
b\ <C 1. As an example, consider the model used in the first studies of so-called "kinetic 
ballooning modes" (in this work, called FLR ballooning modes) with lu ~ u>* p i [? ? ? 
] [normalized as in equation (6)]: 

= <5<J" + lu(uj -uj^S^s-VS^s. (A.6) 

In the limit \s9\ 3> 1, where V oc l/(s9) 2 vanishes, equation (A.6) describes harmonic 
waves. In contrast, the solutions (A. 3) have the following notable properties: 

• The parallel "wavelength" (~ and both the phase and the group velocity 
are all proportional to 9; whereas, in the lowest-order FLR model, they are 
constant. 

• The wave amplitude exhibits algebraic growth proportional to l^l 1 / 2 . The 
exponential increase/decrease in the magnitude of damped/growing solutions (due 
to Im{fc||} oc 7 Sj 0) known from the lowest-order FLR model is also turned into 
algebraic increase/decrease in the full FLR model. 

These differences can be seen in figure 3 which shows shooting code results for aTAE 
trapped in the central potential in (a) the ideal MHD limit [equation (A.6) with 
w* P i = 0] and (b) the full FLR MHD model [equation (7)]. Here, the bound state 
component is located in the region \9\ < tt and the outgoing continuum wave in the 
region |0| > n. 

After Fourier transformation, the mode structure Sip(nq — m) exhibits singular 
spikes at radii satisfying \nq(r) — m\ 2 = u) T (u> T — u)* p i) > 0, where uJZpl is the effective 
diamagnctic frequency defined in equation (19). In FLR MHD, the Alfvcn velocity 
diverges as \9\ — > oo, because the ion polarization current vanishes (and there is no 
electron inertia or collisionality in our model to replace it), so the locations of the 
singularities move to nq{r) — m = 0. The effect is illustrated in figure A2. Note that, 
since the propagation velocity is finite, |0| — *■ oo translates to t — * oo, which means 
that the continuum frequency sweeps up for Re(fcii ) =/= 0. 



Implications for the construction of quadratic forms 

The algebraic increase in the wave amplitude seen in figures 3(b) and 4 is consistent 
with the analytical result (A. 3), so it is not a numerical artifact. However, this raises 
the question whether it is possible to construct physically meaningful quadratic forms 
for a model, which is not physically self-consistent for |0| 3> l/|fcois|- In other words, 
the question is whether and how square-integrability is satisfied when these solutions 
are destabilized by kinetic effects, as is done in [? ]. As will be shown in this section, 
the answer lies in the fact that, unlike in the ideal-MHD limit given by equation (17), 
the magnitude of | <5NI/ S | 2 is not an appropriate measure for the wave energy density. 
Written in quadratic form, equation (7) becomes 

i$ b = 5C; (A.7) 
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nq-m nq-i 



Figure A2. Radial mode structure 5tj){nq — m) ( ) of an aTAE(l, 0) in the 

case with s = 1.0 in table 1 for a = 2.7. The ideal MHD (a) and the FLR 

MHD solution (b) is shown. The mode structure exhibits singularities ( ) 

at the locations where the effective mode frequency ^/o; r (u; r — u> tp i) ( — • — ) 

matches the frequency of the local continuum ( ). In (b), the continuum 

frequency ( ) diverges as m ax — > oo (which translates to t — * oo due to finite 

propagation velocity), but it is plotted here for a finite domain size to fit the 
numerical result for Sip(nq — m). 



7 max) u m&x\ ? 



where, for the simulation domain 8 G [— (9 n 

$ b = i [£v|> s W/_7; M . (A.8) 
dO \\5%f + V\5V B f] (A.9) 



/ dQ (Kn,c + Kn,r + K.FLr) |<$* 8 



Here, $b measures the energy flux through the boundary and 5C is the Lagrangian (see 
equation (31) in [? ]). The structure of V c g and ^ indicates that \Sip\ 2 = \S^ S \ 2 / f 
represents the wave energy density, even though is the relevant wave function; 
i.e., a solution of the SAW Schrodingcr equation (7). 

Our simulations using AWECS [? ] confirm that |<5"0| 2 decays like with x w 1 

for marginally stable modes and x > 1 for unstable modes. Thus, unstable modes are 
square-integrable and it is meaningful to use quadratic forms for theoretical arguments 
and for analysis of simulation results, which is done in a companion paper [? ]. 

The above implies that the appropriate normalization constant for the quadratic 
forms is J d9\S^> s \ 2 / f [as in equation (21)]. Although the normalization constant is 
irrelevant for any single case, the choice is crucial for plotting meaningful results for 
parameter scans [? ]. 

Note that the first line of equation (A.9) corresponds to 28W{ which measures the 
ideal MHD potential energy and the cu 2 term on the second line measures the ideal 
MHD inertia. However, since the eigenvalue of an oTAEs is determined in the region 
between the turning points of the potential barriers, there is no two-scale structure 
and no distinction can be made between an "ideal region" and an "inertial layer" (as 
used, e.g., in [? ]), so we prefer not to separate kinetic and potential energies and keep 
them unified in the Lagrangian $£[?]. 

Implications for numerical simulation of a linear initial-value problem with aTAEs 

Numerical simulations are always done in a finite-size domain, so outward propagating 
waves must be absorbed by appropriate boundary conditions in order to avoid 



a TAE: I. Ideal MHD and FLR effects 



23 



unphysical reflections. Here, we discuss two methods which may be used in linear 
gyrokinctic initial-value simulations of ojTAEs: 

• Matching: In our shooting code, we successfully match the numerical solution to 
a known analytical form given by equation (A. 3) at the boundaries 6 = ±# max of 
the computational domain. This is not as easily done in gyrokinctic initial value 
simulations, since the boundary condition given by equation (A. 3) is a function of 
the eigenvalue w, which is not known a priori. One may attempt to overcome this 
constraint by an iterative or adaptive approach which, however, may adversely 
affect the stability of numerical finite-difference and integration schemes, and 
appears to require additional artificial damping. 

• Absorbing boundary: A popular method used in initial value codes is the absorbing 
boundary condition, where one imposes an artificial damping rate in the outer 
reaches of the simulation domain. However, this is effective only for those 
components of the wave, which have wavelengths smaller than the damping 
region. For the dominant long-wavelength component of the FLR continuum 
waves described by equation (A. 3), such an absorption layer acts like a reflecting 
wall. 

These examples illustrate that it is difficult to treat FLR continuum waves 
accurately in initial- value codes. Fortunately, it is often possible to tolerate unphysical 
reflections off the boundary; in particular, when the following constraints are satisfied: 

• The mode is driven strongly unstable. By the time the continuum wave returns, 
its energy is exponentially small compared to the bound state. 

• The mode is unstable and deeply trapped between large potential barriers, so the 
fraction of returning wave energy which penetrates the barriers is exponentially 
small. 

Clearly, the most challenging case is the simulation of weakly bound states near 
marginal stability. 



